#!/usr/bin/perl -w
use strict;
my $window = 0;
use Bio::DB::SeqFeature::Store;
my $dir = shift || 'gff3';
my $domains = shift || 'batrachochytrium_dendrobatidis_JEL423_1.InterPro.tab';
my $db =  Bio::DB::SeqFeature::Store->new
  (-adaptor => 'berkeleydb',
   -dir     => $dir);

open(my $ifh => $domains) || die $!;
my %gene_domains;
my %go;
while (<$ifh> ) {
  chomp;
  my @row = split(/\t/,$_);
  $gene_domains{$row[0]}->{join(" ","$row[3]:$row[4]",$row[5])}++;
  my $last = pop @row;
  next if $last !~ /GO/;
  $go{$row[0]}->{$_}++ for split(/,\s*/,$last);
}
my @features   = $db->get_features_by_type('gene:Broad');
print join("\t", qw(GENE SNP_COUNT DOMAINS GO)), "\n";
open(my $fh => ">SNP_all.tab");
for my $f ( @features ) {
  my @SNPs = $db->features( -type => ['SNP:Illumina_BWA'],
			    -seq_id => $f->seq_id,
			    -start => $f->start - $window,
			    -end   => $f->end + $window,
			  );
  print join("\t",$f->name, scalar @SNPs,
	     join("; ", keys %{$gene_domains{$f->name} || {}}),
	     join("; ", keys %{$go{$f->name} || {}}),
	    ), "\n";
  for my $snp ( @SNPs ) {
    print $fh join("\t", $f->name, $snp->seq_id, $snp->start,
		   $snp->end, $snp->get_tag_values('RefAllele'), 
		   $snp->get_tag_values('Alleles')),"\n";
  }
}
